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ABSTRACT 

Symplectic integrators are widely used for long-term integration of conservative astrophysical prob¬ 
lems due to their ability to preserve the constants of motion; however, they cannot in general be 
applied in the presence of nonconservative interactions. In this Letter, we develop the “slimplectic” 
integrator, a new type of numerical integrator that shares many of the benefits of traditional sym¬ 
plectic integrators yet is applicable to general nonconservative systems. We utilize a fixed-time-step 
variational integrator formalism applied to the principle of stationary nonconservative action devel¬ 
oped in Galley (2013); Galley, Tsang, & Stein (2014). As a result, the generalized momenta and energy 
(Noether current) evolutions are well-tracked. We discuss several example systems, including damped 
harmonic oscillators, Poynting-Robertson drag, and gravitational radiation reaction, by utilizing our 
new publicly available code to demonstrate the slimplectic integrator algorithm. 

Slimplectic integrators are well-suited for integrations of systems where nonconservative effects play 
an important role in the long-term dynamical evolution. As such they are particularly appropriate 
for cosmological or celestial N-body dynamics problems where nonconservative interactions, e.g. gas 
interactions or dissipative tides, can play an important role. 


1. INTRODUCTION 

Symplectic integrators are a class of mappings that al¬ 
low for numerical integration of conservative dynamical 
systems and which, up to round-off, exactly preserve cer¬ 
tain constants of motion (e.g. the symplectic form). As 
a result the integrations do not suffer from numerical 
“dissipation” which would cause an unphysical drift over 
many dynamical times. Due to these properties, sym¬ 
plectic integrators are widely used in the long-term inte¬ 
gration of many physical systems, particularly in celes¬ 
tial dynamics (Wisdom & Holman 1991; Gladman et al. 
1991; Levison & Duncan 1994; Rein & Tamayo 2015). 

Conservative variational integrators (see e.g. Marsden 
& West 2001) are a subclass of symplectic integrators 
where the mappings are determined by the extremiza- 
tion of a discretized action. 5 The discretized action can 
inherit the symmetries of the full action such that, by 
Noether’s theorem, the discrete equations of motion ex¬ 
actly conserve the symplectic form and the momenta. 
Since discretizing the time coordinate breaks the con¬ 
tinuous time-shift symmetry, fixed-time-step variational 
integrators that preserve the symplectic form and the 
momenta cannot also conserve energy (Ge & Marsden 
1988). However, the energy error tends to be bounded 
by a constant, even over long integration times (see e.g. 
Lew et al. 2004, and references therein), in contrast with 
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traditional integration methods where the error tends to 
grow with time. 

Variational integrators can be applied to some dis¬ 
sipative problems using the Lagrange-d’Alembert ap¬ 
proach (Marsden & West 2001; Lew et al. 2004). Here, 
we utilize the more general nonconservative action princi¬ 
ple, recently developed by Galley (2013); Galley, Tsang, 
& Stein (2014). This formalism was developed to accom¬ 
modate generically the causal dynamics of untracked or 
inaccessible degrees of freedom that might result from an 
integrating-out or coarse-graining procedure at the level 
of an action/Lagrangian/Hamiltonian. 

In this Letter, we develop variational integrators from 
the nonconservative action principle. The resulting map¬ 
pings are no longer symplectic, as the symplectic form 
(and momenta) are no longer conserved, but evolve ac¬ 
cording to the nonconservative dynamics. We instead 
refer this type of numerical integrator as “slimplectic” 
since phase space volumes tend to slim down for dissipa¬ 
tive systems. We show that our method inherits many 
of the same performance features of the symplectic in¬ 
tegrator. Previous works have demonstrated some suc¬ 
cess by including weakly dissipative forces in 2nd-order 
symplectic integrators (see e.g. Malhotra 1994; Cordeiro 
et al. 1996; Mikkola 1997; Hamilton et al. 1999; Zhang 
& Hamilton 2007); in fact, these can be shown to be 
particular cases of our more general (arbitrary order) 
method. For brevity, we focus on a basic fixed-time- 
step slimplectic integrator leaving further developments, 
such as adaptive time-stepping and detailed discussion of 
Noether current evolution, to a longer follow-up paper. 


2. NONCONSERVATIVE LAG RAN GIAN MECHANICS 

Nonconservative Lagrangian mechanics accommodates 
nonconservative interactions and effects by first formally 
doubling the degrees of freedom, q -A (gi, # 2 )- The action 
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describing the dynamics of these doubled variables is 


s = j Hqi,i,qi,i,t)dt ( 1 ) 

where the (nonconservative) Lagrangian is 
A(<? 1 , 2 , qi, 2 , t) = L(q 1 ,qi,t)-L(q 2 , q 2 , t)+K(q h2 , qi, 2 ,t). 


L is the usual Lagrangian, which is an arbitrary function 
of coordinates, velocity, and time and describes the con¬ 
servative sector of the system (i.e., dynamics in the ab¬ 
sence of nonconservative effects). However, K is another 
arbitrary function that couples the variables together, 
vanishes when q\ = < 72 , and accounts for any generic 
nonconservative interaction. Note that A is completely 
specified once L and K are given. 

After all variations of S are performed the two vari¬ 
ables are identified with each other, q\ = q 2 = g, which 
is called the physical limit (PL). In some cases, it is con¬ 
venient to work with more physically motivated coordi¬ 
nates, g_ = qi — q 2 and q+ = (#1 +^)/2. The former can 
often be considered like a virtual displacement and van¬ 
ishes in the PL while the latter is the surviving physically 
relevant combination. 

Requiring S to be stationary under variations of the 
doubled variables and taking the PL leads to nonconser¬ 
vative Euler-Lagrange equations of motion, 


d dA 
dt dq- 


8A 

dq- 


= 0 


PL 


( 2 ) 


or in terms of L and P, 
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Figure 1. The time discretization used for the Galerkin-Gauss- 
Lobatto (GGL) variational integrator. The action is approximated 
using the Gauss-Lobatto quadrature method, where for each in¬ 
terval At, the area under the nonconservative Lagrangian is ap¬ 
proximated by the weighted sum of the Lagrangian, A(q±, q±,t), 
evaluated at each of the (r + 2) quadrature points. This quadra¬ 
ture rule (and the variational integrator) is accurate up to order 
(2 r + 2). The generalized velocities q± are approximated using 
the derivative of the (r + l)th-order interpolating function at each 
quadrature point. 


where = —1, x r +\ = +1, and Xi (for i G {1 .. .r}) is 
the ith root of dP r+ i/<P, the derivative of the (r + l)th 
Legendre polynomial, P r+ \{x). For a given nonconser¬ 
vative Lagrangian functional A(g±, g±, £), we can ap¬ 
proximate the degrees of freedom q n ,±(t) = (j) n ,±{t) + 
0(At r+2 ) using the cardinal-function interpolation for 
this choice of quadrature points. 

We then have the approximation q n ,±(t) ~ 0 n ,±(£)> 
which can be conveniently evaluated at the quadrature 
points using the derivative matrix (see e.g. Boyd 2001, 
2015), 


There are multiple ways to specify P, which depend on 
the particular problem in question. More details about 
this aspect and the nonconservative action formalism in 
general can be found in §11 of Galley, Tsang, & Stein 
(2014), including the evolution of Noether currents ac¬ 
cording to nonconservative processes described by P. 

3. VARIATIONAL INTEGRATORS FOR 
NONCONSERVATIVE SYSTEMS 

Variational integrators numerically approximate the 
behavior of a system by implementing the exact 
equations of motion for a closely-related discrete ac¬ 
tion (Marsden & West 2001; Brown 2006). Here we will 
apply it to the nonconservative action principle described 
above. 

To construct a variational integrator we need to make 
a choice of discretization of the action integral S = 

Jtf A (q±,q±,t)dt. 

A choice that provides a time-reversal symmetric 
discretization (thus an even order method; Farr & 
Bertschinger 2007) is Gauss-Lobatto quadrature, illus¬ 
trated in Figure 1. On the time interval t G [£ n ,£ n+ 1 ] 
with At = £ n+ 1 — t n , we have the set of r + 2 quadrature 

points t n , {tnh r i= i,t n +\, with 

A* = L + (1 + (4) 


Dij — < 


—(r + l)(r + 2)/(2Ai) 
(r + l)(r + 2)/(2At) 
0 

2P r+ i(xt) 

. P r +i(xj){xi - Xj)At 


i = 3 = 0 
i = j = r + 1 
* = 3 + 0, r + 1 


* + 3 


such that 


r+1 


A,±(4^) = F Di i 4n,± = A*± 

j =0 


( 5 ) 


(6) 


where for notational compactness we define tn^ = t n , 
tn +1) = t n + 1 , and q^± = q±{t ( n ) ). 

Using Gauss-Lobatto quadrature, any integral func¬ 
tional f Fdt [for example F G {A, L,P}] has a discrete- 
quadrature approximation on the time interval [t n , t n + 1 ]. 


This discrete functional F£ is 6 


Td(^n,±5 {^n,±}i= 1 J (Zn+l*±j t n ) — ^ ^ w iF(q n ,±i ^n,±J L h 

i =0 


r+1 


(i) t (i)j 


= F, 


d > 


( 7 ) 


For L^ we can drop the d=, indices. 
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where the Gauss-Lobatto quadrature weights Wi are 
given by 


At 


Wi 


(r + l)(r + 2)[P r+1 (xi)} 2 ' 


(8) 


Now we approximate the action over an interval 
[to,tN+ l]j 

rtN + l 

S[t 0 ,t N+1 ] = / A (q±,q±,t)dt 

Jto 

= <Sd[t(htN+i] + (9(At 2r+3 ), (9) 

where the discretized action is defined as 

N 

tjV+l] = ^ ^ A^(ff n; ±, (?n+l,±; ^n)- (10) 

n=0 


We refer to this discretization choice as the Galerkin- 
Gauss-Lobatto (GGL) method (Farr & Bertschinger 
2007). 

The discrete action Sd[to,tN+i\ from (10) can then be 
extremized over values q n5 _ and , and the physical 
limit imposed, to generate the discretized equations of 
motion for each n G [1,7V], 


r 1 

. _ 
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n 
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dq n . 


PL 
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= 0, 

(11a) 

= 0, 

(lib) 


since from Figure 1, we see that each q^_ only con- 
tributes to a single in the discretized action (10), 
while each q n _ appears in both A^ _1 and A^. 

In terms of L^ and the equations of motion are 


9LT 1 | 
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9K2- 1 

dq„- 
dL r ‘ 


dq • 


(i) 


+ 


dK2 

dq n - 

' dK2 ' 

Mu}-. 


= 0 


PL 


= 0 . 


PL 
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We now introduce the discrete momenta 7r n , defining the 
nonconservative (slimplectic) GGL variational integrator 
map (g n ,7Tn) (q n +u tRiTi), by splitting the equation 
of motion (12a) into 
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PL 


Given initial values of (g n ,7r n ), the values of q n +1 are 
determined implicitly by (13a), while the values for the 

q$ intermediate points are given implicitly by equa¬ 
tion (13c) for i G {1 ...r}. The final momenta 7r n+ i 
can then be determined explicitly from (13b). 



Time, t [( m/kf 2 ] 



Figure 2. Top: Numerical solutions for the damped harmonic 
oscillator, described by L = mq 2 /2 — kq 2 /2 , K = —Xq+q-, with 
A = 10 _4 (mfc) 1 / 2 , and fixed time steps At = 0.1 (ra/fc) 1 / 2 . Initial 
conditions were taken to be g(0) = 1 and g(0) = 0. The 2nd order 
RK solution, shown in green, is unstable for these parameters and 
diverges significantly. The 4th-order RK solution (blue-dashed) 
cannot be readily distinguished from the 4th-order slimplectic solu¬ 
tion (solid-orange) in this plot. The 2nd-order slimplectic solution, 
(solid-red), gives nearly the correct amplitude after ~ 10 5 time 
steps, as the energy evolution is accurately followed, though a phase 
shift is evident. Bottom: The fractional energy error relative to 
the energy given by the analytic solution at each time. We see that 
while initially the RK and slimplectic energy errors are compara¬ 
ble at each order, the RK energy errors grow roughly linearly with 
time, while the slimplectic energy error remains bounded. 


Noether’s theorem for conservative actions can be 
shown to generalize to nonconservative systems where 
the corresponding Noether currents evolve in time due 
to a non-zero K (Galley, Tsang, & Stein 2014). One 
can show that for continuous symmetries of the conser¬ 
vative action, which remain after discretization, discrete 
Noether currents will also evolve due to Kd- Thus, trans¬ 
lational or rotational symmetries, for example, will gen¬ 
erate discrete momenta that evolve according to LQ, up 
to round off and bias error (Brouwer 1937; Rein & Spiegel 
2015). Additional error compared to the physical evolu¬ 
tion is only due to the discretization of the action. 

The GGL discretization does not preserve the time- 
shift symmetry preventing energy evolution from being 
precisely tracked. However, the fractional energy error 
tends to be oscillatory and bounded by a resolution and 
order-dependent constant. 7 We will defer more detailed 
discussion of Noether current evolution to a longer fol¬ 
lowup paper in the interests of space. 

The resulting slimplectic maps are accurate up to or¬ 
der 2r + 2. For r = 0, where no intermediate steps are 
used, the quadrature method is the trapezoid rule, and 
the variational integrator is 2nd-order and equivalent to 

7 Fixed-time step variational integrator methods cannot be both 
symplectic-momentum and momentum-energy preserving (Ge &; 
Marsden 1988), however adaptive time-stepping allows symplectic- 
energy-momentum methods to be developed (Kane et al. 1999; 
Preto & Tremaine 1999; Lew et al. 2003). 
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the Stormer-Verlet “leap-frog” integrator (Wendlandt & 
Marsden 1997). 

It is well known that 2nd-order “leap-frog” integra¬ 
tors can be used for dissipative systems, by inserting a 
dissipative “kick” force into the “kick-drift-kick” ansatz, 
resulting in good energy and momentum evolution prop¬ 
erties. Our approach explains why this simple modifica¬ 
tion works in the 2nd-order system, as it is equivalent to 
the lowest order slimplectic GGL method (see also Lew 
et al. 2004, for a similar Lagrange-d’Alembert approach). 
The slimplectic method allows this to be generalized to 
higher orders and general nonconservative systems 8 . 

4. CODE AND EXAMPLES 

We have developed a simple python code, 
slimplectic, that is publicly available 9 and gen¬ 
erates the fixed-time-step slimplectic GGL integrators 
described above, for use in characterizing the numerical 
technique. The code generates slimplectic solvers of 
arbitrary order (2r + 2) given sympy (SymPy Dev Team 
2014) expressions for L(q,q,t) and K(q±,q±,t). This 
demonstration code is designed to work for arbitrary 
L and IT, and thus has not been optimized as would 
be appropriate for specific problems. In particular, 
the equations of motion (13) are solved with standard 
root-finders, rather than a problem-specific iteration 
method, to be more generally applicable. 

As a basic example in Figure 2 we compare Runge- 
Kutta (RK) and slimplectic integration of a simple 
damped harmonic oscillator, for both 2nd and 4th-order 
methods. Below we also present two basic astrophysi- 
cal examples of non-conservative interactions. All exam¬ 
ples are available as ipython notebooks in our public 
repository. 9 


4.1. Poynting-Roberts on Drag 

We first examine the orbital motion of a dust parti¬ 
cle experiencing Poynting-Robertson drag (Burns et al. 
1979) due to radiation from a solar type star, starting 
with a semi-major axis of 1AU in an eccentric (e = 0.2) 
orbit. The Lagrangian for this system is 


r 1 .9 m GMqTTI 

L=-mq +( 1 -/ ? )^p. 


(14) 


where m is the dust particle’s mass, and q its position. 
The dimensionless factor 


/? = 


3 Lq 

Si xcp GMq(1 


0.058 


_P _ 

cm -3 



(15) 


is the ratio between forces due to radiation pressure and 
gravity, where L 0 and M 0 are the solar luminosity and 
mass, p and d are the density and size of the dust grain. 
The nonconservative potential which generates the cor¬ 
rect Poynting-Robertson drag force is (for |q| « c) found 


8 Similar results can be obtained for Wisdom-Holman-type map¬ 
pings by splitting the nonconservative action into integrable and 
perturbative terms (see e.g. Farr 2009). 

9 The repository for slimplectic is available at http://github. 
com/davtsang/slimplectic. 
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Figure 3. Top: Evolution of the Cartesian x-coordinate for or¬ 
bital motion of a particle experiencing Poynting-Robertson drag 
due to radiation from a solar-type star, with particle density p = 2 
g cm -3 and particle size d = 5 X 10 -2 cm. The particle has ini¬ 
tial semi-major axis of 1AU, and initial eccentricity of 0.2. With 
a fixed-time-step of At = 0.01 yr, the 2nd-order RK integrator 
(green-dashed) is unstable, while the 2nd-order slimplectic integra¬ 
tor (red) behaves much more accurately, despite significant (numer¬ 
ical) precession over a ~1000 year timescale. The 4th-order integra¬ 
tors (RK, blue-dashed; slimplectic, orange) have similar long-term 
amplitude evolution without visible numerical precession. The 
phase errors (not shown) of the RK integrators grow oc £ 2 , while the 
slimplectic integrators have phase error oc t (see Preto & Tremaine 
1999). Bottom: Fractional energy error (compared to a 6th-order 
slimplectic integration) for the system described above. The RK 
errors grow roughly linearly in time, while the slimplectic errors 
are bounded. There is a slight turn-up in the 4th-order slimplectic 
error envelope, most likely due to the build up of round-off, bias, 
or action-discretization error. 

as the virtual work from the known force, 

q+ -q- + A(q+ -q+)(q+ -q-) • 

< 1 + 

(16) 

Methods to determine or derive K are discussed in (Gal¬ 
ley et al. 2014). The first term in square brackets above 
gives the usual drag term, while the second term is due 
to the Doppler shift caused by radial motion. 

The system was integrated for 6000 years using 2nd 
and 4th-order RK (green and blue dashed) and slimplec¬ 
tic (red and orange solid) methods with time-steps of 
At = 0.01 yr. The results and discussion are shown in 
Figure 3. 


K = 


fiGM^m 

2 

cql 


4.2. Gravitational Radiation Reaction 


We also consider two 1.4M 0 neutron stars inspiraling 
from gravitational wave emission. This example demon¬ 
strates (fixed-time-step) integrators for systems where 
the orbital dynamics can change quickly due to noncon¬ 
servative effects. In the post-Newtonian (PN) approxi¬ 
mation, leading-order conservative dynamics for the or¬ 
bital separation, q, are described by Newtonian gravity 


L 





(17) 
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Figure 4. Top: PN radiation reaction, including only the 2.5PN 
dissipative terms through K. Each of the methods (RK and 
slimplectic GGL) are 4th order with fixed timesteps of At = 1000 M 
and At = Tisco- The integration methods blow up when the or¬ 
bital time-scale is comparable to the time step, though the slimplec¬ 
tic method performs significantly better than the RK method for 
the same time steps; for the 1000M time-step the RK method is 
immediately unstable. 

Bottom: Relative error in the orbital radius compared to the an¬ 
alytic adiabatic-approximation solution. 


where M = mi + m 2 = 2.8M Q is the total mass, /1 = 
— 0.7 Mq is the reduced mass, and G = c= 1. 
Dissipative effects from radiation reaction first appear 
at PN order (|q|/c) 5 (or 2.5PN) and are described by K , 
which has been calculated in Galley & Leibovich (2012); 
Galley & Tiglio (2009). After order-reduction, 
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(18) 


where v — fi/M = 1/4. A physically consistent simula¬ 
tion should go to the same PN order in L and K. Here, 
instead, we use the leading order in each as a toy model 
to focus on the numerical method. 

In Figures 4 and 5 we examine the 4th-order RK and 
slimplectic integrators for different choices of time steps, 
At = 1000M and At = Tisco — 92M, the orbital period 
at the innermost stable circular orbit. The initial orbital 
separation is ro = 100M « 414km and corresponds to an 
orbital frequency of ~ 11.5Hz. We compare our numeri¬ 
cal results to analytic solutions to our toy PN equations 



Figure 5. Top: The phase evolution of the integrators for the 
PN radiation reaction depicted in Figure 4. 

Bottom: Absolute orbital-phase errors. In this example, both 
slimplectic integrators (phase error oc £), track the orbital phase 
much better than the equivalent RK integrators (phase error oat 2 ). 


in the adiabatic regime, 



r o 



■r-rffW 

32z/M 5 / 2 


where the orbital period is assumed to be much smaller 
than the radiation reaction time scale for the inspi¬ 
ral (Blanchet 2014). 

All integrators begin to fail when the orbital period 
becomes roughly comparable to the time-step, although 
the slimplectic integrators can get significantly closer to 
this limit than the RK method (the RK integrator with 
time-step At = 1000M immediately becomes unstable). 
In our followup paper, we will demonstrate an adaptive 
time-stepping scheme that will allow efficient slimplectic 
integrations that also precisely evolve the energy. 


5. DISCUSSION 

We have developed a new method of numerical integra¬ 
tion that combines the nonconservative action principle 
of Galley (2013); Galley et al. (2014) with the variational- 
integrator approach of Marsden & West (2001). These 
“slimplectic” integrators allow nonconservative effects to 
be included in the numerical evolution, while still pos¬ 
sessing the major benefits of normally conservative sym- 
plectic integrators, particularly the accurate long-term 
evolution of momenta and energy. 

The discrete equations of motion are found by vary¬ 
ing a discretized nonconservative action and implicitly 
define the slimplectic mapping (g n ,7r n ) -A (g n+ i,7r n+ i). 
Different choices of discretization generate different vari¬ 
ational integrators. Here we have focused on imple¬ 
menting the Galerkin-Gauss-Lobatto (GGL) discretiza¬ 
tion, and demonstrating its long-term accuracy using the 
damped harmonic oscillator, Poynting-Robertson drag 
on a small particle, and a gravitational radiation-reaction 
toy problem. Our results also explain why the modifica¬ 
tion of the 2nd-order “kick-drift-kick” ansatz to include 
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dissipative forces performs so accurately, as this is equiv¬ 
alent to the lowest order version of the slimplectic GGL 
method. 

We have developed a demonstration python code, 
slimplectic , 9 which generates slimplectic integrators 
for arbitrary Lagrangians and nonconservative poten¬ 
tials. Readers are encouraged to test different physical 
systems of interest using this publicly available code, but 
to separately implement problem specific optimizations, 
particularly when solving the implicit equations of mo¬ 
tion. 
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